# required packages - install first
require(estimatr)
require(tidyverse)
require(rdrobust)



###################################################################
###################################################################
#
#    When the state becomes complicit: mayors, criminal actors, 
#   and the deliberate weakening of the local state in Colombia 
#
#                    Camilo Nieto-Matiz
#                 camilo.nieto-matiz@utsa.edu
###################################################################
###################################################################

 

# load dataframe
load("crimcollusion_colombia_main.rda")



################################################################
#                       Tables  
################################################################

#================================================
#                   Table 1
#                Main results  
#================================================

summary(rdrobust(cc_data$tx_0811,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(cc_data$tx_0811,cc_data$share1,c=0,p=2,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))


#================================================
#                   Table 2
#          Effects on political competition 
#================================================

summary(rdrobust(cc_data$comp2011,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(cc_data$parties2011,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(cc_data$cand11,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(cc_data$reelected,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))



#======================================================================
#                              Table 3
#     Mechanisms: eroding property rights and judicial institutions 
#======================================================================

summary(rdrobust(cc_data$cad_0809,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(cc_data$val_0809,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(rdrobust(log(cc_data$inf_0809+1),cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

summary(tx0<-rdrobust(cc_data$just_0809,cc_data$share1,c=0,p=1,  
                      all = T, vce = "nn", kernel="tri",covs=NULL))



#======================================================= 
#                    Table 4
#     State weakening as an non-intentional outcome
#======================================================== 

summary(rdrobust(log(cc_data$para_0811+1),cc_data$share1,c=0,p=1,  
                       all = T, vce = "nn", kernel="tri",  
                       bwselect = "mserd"))

summary(rdrobust(log(cc_data$guer_0811+1),cc_data$share1,c=0,p=1,  
                       all = T, vce = "nn", kernel="tri",  
                       bwselect = "mserd"))

summary(rdrobust(cc_data$royal_0811,cc_data$share1,c=0,p=1,  
                       all = T, vce = "nn", kernel="tri", 
                       bwselect = "mserd"))









################################################################
#                       Figures  
################################################################

#================================================
#                   Figure 1-a
#          Regression discontinuity plot
#================================================

rd.data1 <- cc_data  
# make the binned averages - bins on the left
count <- 1
bin.size <- .02
binx.left <- vector(length=length(seq(-.2, 0, bin.size)))
biny.left <- vector(length=length(binx.left))
last <- -.02
for(j in seq(-2, 0, bin.size)) {
  biny.left[count] <- mean(rd.data1$tx_0811[rd.data1$share1 >= j-bin.size & rd.data1$share1 < j], na.rm=T)
  binx.left[count] <- (j+last)/2
  last <- j
  count <- count + 1
}

# bins on the tight
count <- 1
bin.size <- .02
binx.right <- vector(length=length(seq(0,.2, bin.size)))
biny.right <- vector(length=length(binx.right))
last <- 0.2
for(j in seq(0,2, bin.size)) {
  biny.right[count] <- mean(rd.data1$tx_0811[rd.data1$share1 >= j-bin.size & rd.data1$share1 < j], na.rm=T)
  binx.right[count] <- (j+last)/2
  last <- j
  count <- count + 1
}


temp1 <- data.frame(cbind(binx.left, biny.left))
temp2 <- data.frame(cbind(binx.right, biny.right))
colnames(temp2) <- c("binx.left", "biny.left")
temp <- rbind(temp1, temp2)

# select variables
rd.data1<- select(rd.data1, share1, tx_0811) %>% dplyr::rename(binx.left=share1, biny.left=tx_0811)

# create plot
ggplot(data=temp, aes(x=binx.left, y =  biny.left)) +
  geom_point(data=subset(temp, binx.left >= 0&binx.left <= 2),color="royalblue2",alpha = .4, size = 4.5) +
  geom_point(data=subset(temp, binx.left <= 0&binx.left >= -2),color="tomato2",alpha = .4, size = 4.5) +
  geom_vline(xintercept = 0) +
  geom_smooth(data=subset(rd.data1, binx.left >= 0&binx.left <= .5),method = "loess",formula = y ~ poly(x, 2), se = TRUE, size = 1.5, fill="blue", color = "gray10", span = 1, alpha = .4)  +
  geom_smooth(data=subset(rd.data1, binx.left <= 0.0&binx.left >= -2),method = "loess",formula = y ~ poly(x, 2),se = TRUE, size = 1.5, fill="red", color = "gray10", span = 1, alpha = .4)  +
  coord_cartesian(xlim = c(-.2, .2))+ylim(0,80)+
  labs(x="Margin of Vote", y="Property Taxes") +
  theme_bw()+
  theme(text = element_text(hjust = 0.5, family="Palatino", size=12))
ggsave("tax_rdplot.png", width = 6, height = 4)



#================================================
#                   Figure 1-b
#       RD estimates - different bandwidths
#================================================


# matrix to store homicide estimates
store_bw1 <- as.data.frame(matrix(NA,34,3))
 
# Initiate loop for taxes 08-11
bandw <- sort(c(seq(.04,.2,by=.005), 0.078)) #calonico's optimal bw: 
for (i in 1:length(bandw)) {
  # Estimate model
  summary(bw_loop1 <- rdrobust(y=(cc_data$tx_0811), x = cc_data$share1, c =0, p=1,h=bandw[i],vce = "nn", 
                               kernel="tri", covs=NULL, all = T)) 
  # Extract 3rd row of rd_out() LC function (robust RDD estimate)
  # Place into loop holder
  store_bw1[[i,1]] <- bw_loop1$coef[3]
  store_bw1[[i,2]] <- bw_loop1$se[3]
  store_bw1[[i,3]] <- bw_loop1$pv[3]
  store_bw1$c <- bandw
  store_bw1$c<-round(store_bw1$c, digits = 4)
}
colnames(store_bw1)<-c('coef','se','pv','c')
store_bw1$bw<-ifelse(store_bw1$c==.078,1,0)


# bandwidth plot
store_bw1 %>%
  ggplot()+ 
  geom_errorbar(aes(x = factor(c), ymin = coef - se*1.96,
                    ymax = coef + se*1.96, color=factor(bw), lwd = factor(bw)), width=.4)+
  scale_color_manual(values=c("steelblue3", "blue"))+
  scale_size_manual(values=c(.4,.7))+
  geom_point(aes(x = factor(c), y = coef), lwd = 1.5, shape = 20)+
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)+
  ylim(-90,50)+xlab('') + ylab('') + theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),legend.position = "none")+
  scale_x_discrete(breaks=seq(0.04, .2, by = .02))
ggsave("bw_placebo1_para_tx.png", width = 6, height = 4)







#================================================
#                    Figure 2
#           Effects over time 2012-15
#================================================

# estimate models
summary(t1<-rdrobust(cc_data$tx.2012,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))
summary(t2<-rdrobust(cc_data$tx.2013,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))
summary(t3<-rdrobust(cc_data$tx.2014,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))
summary(t4<-rdrobust(cc_data$tx.2015,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))
summary(t5<-rdrobust(cc_data$tx_1215,cc_data$share1,c=0,p=1,  
                 all = T, vce = "nn", kernel="tri",covs=NULL))

                                                           
# extract coefs and se                                                          
tx_fr <- data.frame(est = c(t1$coef[3],t2$coef[3],t3$coef[3],t4$coef[3],t5$coef[3]),
                    se = c(t1$se[3],t2$se[3],t3$se[3],t4$se[3],t5$se[3]),
                    year = c('2012', '2013', '2014', '2015', '2012-15'))

# plot
tx_fr %>%
  mutate(year = fct_reorder(as.factor(year), as.numeric(as.character(year)))) %>%
  ggplot() +
  geom_errorbar(mapping=aes(x = factor(year), ymin = est - se*1.96,
                            ymax = est + se*1.96), width=.09, color="blue") +
  geom_point(aes(x = factor(year), y = est), 
             lwd = 2, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour = "red", lty = 2, lwd=0.3)+
  ylim(-130,80)+
  ylab('')+xlab('')+theme_bw()+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(hjust = 0),
        axis.text.y = element_text(hjust = 0),
        legend.position="bottom")
ggsave("overtime_tx1215.png", width = 5, height = 4)









#================================================
#                    Figure 3
#       Property taxation in Tarazá and Segovia
#================================================

load("taxes_segovia_taraza.rda")

ggplot(seg_tar, aes(factor(ano),tx, group=municipio, color=factor(municipio),
             size=factor(municipio))) +
  geom_line() + theme_bw() + ylab('Property taxes')+
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        axis.title.x=element_blank()) +
  scale_color_manual(values=c("tomato", "royalblue"))+
  scale_size_manual( values = c(.5,1.5))
ggsave("taraza_segovia.png", width = 5, height = 4)
